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ABSTRACT 

We present a new method for the analysis of pulsar timing data and the estimation of the 
spectral properties of an isotropic gravitational wave background (GWB). By sampling from 
the joint probability density of the power spectrum coefficients for the individual pulsars and 
the GWB signal realisation we can eliminate the most costly aspects of computation normally 
associated with this type of data analysis. We use a 'Guided Hamiltonian Sampler' to effi- 
ciently sample from this higher dimensional (~ 400) space, and show by taking this approach 
we need make no assumptions about the properties of the power spectrum of the GWB, thus 
providing a much more general approach to the problem of pulsar data analysis. When applied 
to the IPTA Mock Data Challenge datasets this allows us to make inferences not only on the 
global properties of the GWB, but also about the individual pulsars in the dataset, whilst also 
providing speedups of approximately two orders of magnitude. 
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1 INTRODUCTION 

Millisecond pulsars (MSPs) have for some time been known to 
exhibit exceptional rotational stability, with decade long observa- 
tions providing timing measurements with accuracies on a par with 
atomic clocks (e.g. |Kaspi, Taylorr& Ryba 1994 Matsaki s, Taylor,| 
|& Eubanks|1997[ >. Such stability lends itself well to the pursuit of a 
wide range of scientific goals, e.g. observations of the pulsar PSR 
B1913-I-16 showed a loss of energy at a rate consistent with that 
predicted for gravitational waves ( Taylor & Weisberg 1989), whilst 
the double pulsar system PSR J0737-3039A/B has provided precise 
measurements of several 'post Keplerian' parameters allowing for 
additional stringent tests of general relativity ^Kramer et al.|2006} . 

By measuring the arrival times (TOAs) of the radio pulses to 
high precision it is possible to construct a timing model: a determin- 
istic model that describes the physical properties of the pulsar e.g. 
its binary period and spin evolution, its trajectory, post-Keplerian 
terms and so on. A detailed description of this process is available 
in the Tempo2 series of papers (Hobbs, Edwards, & Manchester 
[20061 [Edwards, Hobbs, & Manchester 2006 Hobbs et al.,2009j . 
The timing model can then be subtracted from the TOAs result- 
ing in a set of residuals which contain within them any physics not 
correctly accounted for by the timing model. 

In this paper we will be concerned with extracting informa- 
tion from these residuals that results from time-correlated stochas- 
tic signals. These can include additional red noise terms due to ro- 
tational irregularities in the neutron star ( Shan non & Cordes|2010^ 
or correlated noise between the pulsars due to a stochastic gravi- 
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tational wave background (GWB) generated by, for example, co- 
alescing black holes (e.g. Jaffe & Backer"2003' Phinney 2001 1 
or cosmic strings (e.g . Kawasaki , Miyamoto, & Nakayama 2010 
Olmez, Mandic, & Siemens 2010f which could be detected using 
a pulsar timing array (PTA), a collection of Galactic millisecond 
pulsars from which the cross correlated signal induced by a GWB 
could be extracted. Current methods for the analysis of PTA data 
are for the most part extremely computationally expensive. This 
is particularly true for existing Bayesian methods jvan Haasteren[ 
[erar] ( [2009l l, |van Haasteren & Levin[pOT2t henceforth vH2009, 
vHL2012) with large dense matrix inversions resulting in a scaling 
factor of approximately O(w^). Recently a newer method has been 
proposed in [van Haasteren[j2012i l where lossy data compression is 
used to reduce the length of time these matrix inversions require, 
resulting in a speed up of ~ 3-6 orders of magnitude over pre- 
vious methods. This method however results in extreme memory 
usage, and was only applied in the case where two free parameters 
were fitted (i.e. the amplitude and slope of a power law fitted to the 
GWB). 

In this paper we present an alternative approach to perform- 
ing a Bayesian analysis on PTA data that results in a speed up of 
approximately two orders of magnitude, and scales as 0(nlog(n)) 
that is not limited by the number of free parameters fitted, and uses 
~ 1GB of system memory for the analysis of the IPTA Datasets. 
We accomplish this by rephrasing the likelihood to eliminate all 
matrix-matrix multiplications and costly matrix inversions, replac- 
ing them with matrix-vector operations and sparse, banded matrix 
inversions, whilst still retaining the ability to make robust statis- 
tical inferences about the white and red noise present in the PTA 
data with the same precision as in vH2009/vHL2012. We eliminate 
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the need for specifying any prior form for the shape of the corre- 
lated power spectrum induced by a GWB, or the red noise present 
in a particular pulsar at the point of sampling, and provide a means 
for easily detecting bright single source objects that might affect 
any method that attempts to fit a prescribed model for the power 
spectrum, and individual pulsars that are statistical outliers within 
the PTA caused by unknown local effects. We perform the sam- 
pling process using a Guided Hamiltonian Sampler (GHS) (Balan, 
Ashdown & Hobson in prep, henceforth B12) which provides an 
efficient means of sampling over large numbers of dimensions (po- 
tentially > 10"^). 

In Section[2]we describe the derivation of the new likelihood, 
in Section|3]we describe the sampling method used, and in Section 
|4]we describe the results of this method as applied to the first and 
second datasets from the Open IPTA Data Challenge, and make 
comparisons with the methods described in vHL2012. 



2 ESTIMATING THE POWER SPECTRUM 

For any pulsar we can write the TOAs for the pulses as a sum of 
both a deterministic and a stochastic component: 



(1) 



where t^„^ represents the n TOAs for a single pulsar, with ijet and 
<sto the deterministic and stochastic contributions to the total re- 
spectively, where any contributions to the latter will be modelled 
as random Gaussian processes. In estimating the timing model pa- 
rameters for the pulsar, a standard weighted least-squares fit as per- 
formed in packages such as Tempo2 will model the stochastic con- 
tributions purely as white noise characterised by the TOA uncer- 
tainties. Thus any contribution to /jto not described by the TOA 
uncertainties will be absorbed by the timing model fit and so when 
the timing model is subtracted from the data, any attempt to char- 
acterise the power spectrum of the resulting residuals will be incor- 
rect. 

In order to account for this, we begin by following the ap- 
proach of vHL2012 which we describe in brief here so as to aid 
subsequent discussion. A set of pre-fit timing residuals St^^^ are 
produced using an initial estimate of the m timing model parame- 
ters ySo, . From here a linear approximation of the timing model can 
be used such that any deviations from the initial guess of the timing 
model parameters are encapsulated using the m parameters 6, such 
that 



(2) 



This allows us to write the timing residuals St in this linear approx- 
imation to the timing model as: 



6t = Stp,^ + Me 



(3) 



where M is the n x m 'design matrix' which describes the depen- 
dance of the timing residuals on the model parameters. 

We can then write the Bayesian likelihood for the timing resid- 
uals as (vH2009): 



P{St\e,<l>). 



^ :exp(--(dt - MefC'\6t - Me)] (4) 



V(2;r)"detC \ 2 



where the n x n covariance matrix C describes the stochastic con- 
tributions to the timing residuals such that 



and is described by a set of parameters ^. 

We can then marginalise over all variables e in order to cal- 
culate the probability for a particular set of parameters for the 
stochastic contributions to the residuals, i.e. 



I'" 



^ exp\--(6t - MefC-\6t - Me)|.(6) 
V(2;r)'"detC \ 2 / 



In vHL2012 this marginalisation is performed analytically to 



give: 



-^(2;r)(''-"')det(G^CG) 



exp i-l-6fG(G^CGr^G^6t\ (7) 



where G is an positive definite symmetric nx (n - m) matrix, the 
derivation of which will not be described here. 

For the IPTA Data Challenge, data sets consisted of 130 resid- 
uals for 36 pulsars such that n = 4680. G therefore is ~ 4500x4500, 
and so the bottleneck in this calculation comes from the matrix 
inversion that must occur for every likelihood calculation, along 
with the set of matrix-matrix multiplications required to calculate 
G'^CG. 

Our goal is to remove this obstacle by rephrasing the likeli- 
hood such that 1) no matrix-matrix multiplications and 2) no com- 
putationally intensive (i.e. 0(«^)) matrix inversions are required in 
the evaluation whilst retaining the ability to accurately determine 
the power spectrum of the stochastic contributions to the residuals 
in these datasets, with a scaling of 0(;2log(n)). 

We do this first by writing our timing residuals dt as the sum 
of a signal s that we are interested in parameterising, and some 
additional white noise n so that we have 



6t = s + n. 



(8) 



We can expand i in terms of its Fourier coefficients a so that 
5 = Fa where F denotes a Fourier transform. For a single pulsar the 
covariance matrix (f of the fourier coefficients a will be diagonal, 
with components 



■■ (fiSij, 



(9) 



where there is no sum over /, and the set of coefficients {(fj] repre- 
sent the theoretical power spectrum for the residuals. When dealing 
with a signal from a stochastic gravitational wave background how- 
ever it is crucial to include the cross correlated signal between the 
pulsars on the sky. We do this by using the Hellings-Downs relation 
l |Hellings & Downs|1983"l >: 

3 1-cos(6',„„) /I - cos(e,„„)\ ll-cos(0,„„) 1 1 
= 2 2 2 n 2 + 2"'2'''""^^°^ 

where 9,„„ is the angle between the pulsars m and n on the sky and 
a,,,,, represents the expected correlation between the TOAs given an 
isotropic background. With this addition our covariance matrix for 
the fourier coefficients becomes 



■ Ct,„„<pi6ij 



(11) 



where there is no sum over which results in a band diagonal ma- 
trix for which calculating the inverse is extremely computationally 
efficient. 

We then sample from the joint probability density of the power 
spectrum coefficients and the signal realisation Pr(j^,j,a | 6t), 
where here a refers to the concatenated vector of all coefficients 
cij for all pulsars, which we can write: 



(5) 



Pr({,^,l,fl|*0 oc Pr(<5<|fl) Pr(a|{<^,l) Pr((^,l) 



(12) 
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and then marginalise over all a in order to find the posterior for 
the parameters {(fj] alone. For our choice of Pr(j^, j) we use a uni- 
form prior in log^ space as the scale of the coefficients is largely 
unknown below some upper limit, and draw from the parameter 
Pi = log((/),) instead of i/?, which has the added advantage that we 
avoid unnecessary rejections due to samples which have negative 
coefficients in the sampling process. Given this choice of prior the 
conditional distributions that make up Eq.[T2]can be written: 



Pi(St\a) oc 



1 



^det(G^NG) 



exp 



{6t - FaYGiG^mr^G^iSt - Fa) 



(13) 



where N = (nn^^ and represents the white noise errors in the resid- 
uals, which follows from Eq.|7]with N replacing C, and substituting 
St - Fa for 6t, and: 



Pr(a|jp,l) oc 



1 



: exp 



*T -1 

a if a 



(14) 



Note that we can calculate G(G^NG)"'G'^ before the sampling 
starts and store it in memory which eliminates the need for any 
dense matrix inversions, or matrix multiplications within the likeli- 
hood calculation. It is important to state that this does not mean we 
are unable to fit for any of the properties of the white noise. Using 
the fact that N is a diagonal matrix, and G are block diagonal and 
obey the relation: 



(15) 



the inverse term immediately simplifies to N'"' . If we wish to model 
the noise for each pulsar independently, where each pulsar p has Op 
residuals, nip model fit parameters and cr,, white noise amplitude, 
N' will be a diagonal matrix where the first (oi - nii) entries are 
equal to cTj, the next (02 - m^) entries will be equal to cr^ and so 
on. Again exploiting the block diagonal nature of G we can then 
rearrange: 

GN' 'G^ = N 'GG'^ (16) 

where we now have our original nx n matrix N where the first 0| 
entries are equal to a\, the next 02 entries will be equal to a\ and 
so on. For a total of np pulsars in the PTA data we can therefore 
rewrite Eq|I3|as: 



vm\aA<Tp}) 



n 



exp 



1 



I 



2 " 



{6tp - FpapY Gp 



Gl(Stp - Fpap)] , 



(17) 



where the subscript p refers to an individual pulsar in all cases. 
Note that we are still calculating the determinant for N',,, which is 
simply CTp'' rather than N,,, and thus we need only store GG^ in 
memory and we have eliminated any matrix-matrix multiplications, 
whilst still retaining the ability to fit for the noise on a per pulsar 
basis. As with the power coefficients <p we instead draw from the 
log uniform distribution cr,, = 10^^'' and sample uniformly in log 
space. 

For the sake of simplifying our notation we therefore redefine 



N 



GN' 'G' 



(18) 



We then draw the samples from the joint space described in 
equation Eq[T2] using a GHS (B12) which is described in the fol- 
lowing section. 



3 GUIDED HAMILTONIAN SAMPLING 

For a detailed account of both Hamiltonian Monte Carlo (HMC) 
and GHS refer to BI2, here we will describe each only in brief. 
(HMC) sampling ( [Duane et al.|I987l l has been widely applied in 
Bayesian computation (Neal 1993), and has been successfully ap- 
plied to proble ms with extremely large num bers of dimensions 
(~ 10* see e.g. [Taylor, Ashdown, & Hobson| l|2008 1). We define 
a 'potential energy' which is related to our posterior distribution 
Pr(jc) by: 



T(jt:) = - ln(Pr(jt:)) 



(19) 



where x is the A' dimensional vector of parameters to be sampled. 
Each parameter x; must be assigned a mass m, and a momentum p, 
so that we can write our Hamiltonian as: 



(20) 



The sampler is given a start point (x, p) and a set of momenta 
which are drawn from a set of A' uncorrelated Gaussian distribu- 
tions of width nii in dimension The system can then evolve deter- 
ministically from then for some length of time t using Hamilton's 
equations 

dx, _ dH 

At dpi ' 

dpi^ _ _dH_ _ d'Vix) 

dxj 



df 



dx, 



(21) 



(22) 



After it has reached its new position {x' ,p') that point will be 
accepted with a probability 



p = min [1, exp(-(5//)] 



(23) 



where SH - H(x',p') - H(x,p). A new set of momenta can then 
be drawn and the process repeats. 

In order to perform the integration along the systems trajectory 
at each state we use a 'leapfrog' method as is common practice. 
Here steps are taken of size A such that n^A = t such that: 

A d'¥(x) 



P,(t) - 



2 dxj 



Xiit + X) = Xiit) + 



Pi (t + X)' 



A\ Ad'Vix) 
2 dxi 



(24) 



(25) 



(26) 



i(r+.l) 



until f = T where r is varied to avoid resonant trajectories. HMC 
thus requires a large number of adjustable parameters, the mass m,, 
step size A, and the number of steps in the trajectory. Adjusting 
the step size or the mass produces similar effects (Neal 1996) and 
so one is usually fixed and the other tuned during sampling. 

GHS eliminates much of this tuning aspect by using the Hes- 
sian H of the joint probability distribution calculated at its peak to 
set the step size A for each parameter. The masses m, are then set to 
unity and the only tuneable parameter that remains is a scaling pa- 
rameter rj which is chosen such that the acceptance rate for the GHS 
is ~ 68% (see B 12). Once the Hessian has been calculated, one then 
determines its N eigenvalues w, and N normalised eigenvectors e, . 
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Denoting the matrix containing tliese normalised eigenvectors as 
its columns by S, one first defines a new set of variables x' = S^x 
in which the Hessian becomes diagonal with the eigenvalues ai, as 
its diagonal entries. One then rescales each x'. to obtain a new set of 
variables y, = ^^fjjj x'./rj. It is straightforward to show that the new 
variables are related to the original variables by 



flmax = (F^N"'F + ^-'r'F^r'5f, 



(38) 



(27) 



Consequently, in the new variables, the Hessian at the peak has the 
trivial form rf\. One then performs Hamiltonian sampling employ- 
ing the standard leapfrog method, but in terms of the new variables 
y, rather than x. Thus, GHS may be considered simply as standard 
HMC, but performed in a set of variables (or coordinates) that are 
tailored to the target distribution, namely the scaled eigendirections 
of the Hessian at its peak. Consequently, although GHS can take ad- 
vantage of a scenario where Pr(jt:) possesses a single well-defined 
peak (with zero gradient), it does not rely on this, since it retains 
the generality of standard HMC. 

Rather than working in terms of the new variables y, one can, 
if desired, return to using the original variables x, in which case the 
relation l[27j shows that the leapfrog steps take the modified form 



phi) 



Pit) - -T/H-'/^V,^ 



x(t + A) = x{t) + ArjH 



-1/2 



(28) 



(29) 



(30) 



Note that whilst performing the eigen-decomposition can be 
costly, it is only necessary to do so once at the start of the sampling 
where it takes ~ 10s and is thus negligible within the scope of 
the sampling process as a whole. Therefore in order to perform 
sampling we need the following: 

• The gradient of for each parameter x, 

• The peak of the joint distribution 

• The Hessian at that peak 

The gradients of our parameters are given by the following: 
da 

1 / -i'5v\ 1 ,T -,dip 
^r-l'T %)"2« ^ 
and the components of the Hessian are: 

a^ 
a^ 
asf 

a^ 



ln(10)(o, - m,) - \niW){5ti - FiUifN \sti - FiUi) 



F^N F + (p- 



21n(10)^((5 /, - FiUifN \dti - Ftad 



.T -1 9^ -I 9ip _i n « -r -1^ -1 
-a V -^V -^V a - 0.5a —^p 'a 
dp, dpi dp- 



a^f _ _ _id<f_ 
dpida ^ dpi^ 



(31) 
(32) 
(33) 

(34) 
(35) 

(36) 

(37) 



For a set of power spectrum coefficients {p,) and white noise 
coefficients jE, j we can solve for the maximum set of Fourier coef- 
ficients a^ax analytically. Rearranging Eq.[3T|we find: 



so when searching for the global maximum we need only search 
over the subset of parameters {pj. This is achieved by using a parti- 
cle swarm algorithm (Kennedy & Eberhart 1995;2001 and for uses 
in cosmological parameter estimation see e.g. [Prasad & Souradeep| 
( |2012 t, and for a description of the particle swarm method applied 
to PTA data in this context see (Taylor, Gair & Lentati 2012) to 
efficiently find the maximum in ~ 1 minute using ~ 5-10 cores 
per free parameter. We now apply this method to the IPTA Open 
Datasets 1 and 2. 



4 RESULTS 

4.1 Open Dataset 1 

The IPTA Data Challenge 1 first open dataset consists of a set of 
130 timing residuals for 36 pulsars and contains only white noise 
with amplitude cr = 100ns for each pulsar, and an injected GWB 
power spectrum with a characteristic strain spectrum given by: 



hcif) = A„ 



f 



lyr- 



(39) 



with Al, a dimensionless amplitude at a frequency of (yr ') and a 
a power law index. The parameters for the injected spectrum are 
A,, = 5 X lO-'"* and a = -2/3. 

Parameterising the spectral density as in vHL2012: 



S(f) 



1 



lyr ' / \ lyr 



/ 



(40) 



the strain spectrum will result in an observed spectral density within 
the residuals of: 



S(f) 



12;r2 ' \lyr-i 



-13/3 



(41) 



Defining f„ = n/T where T is the total observation length of 
the dataset, we fit for the Fourier coefficients corresponding to some 
set [n]. The power spectrum for each of the Np pulsars is therefore 
described by 2« -I- 1 free parameters; the vector a containing the 2n 
Fourier coefficients (n real, n imaginary) and the parameter E,, de- 
scribing the log amplitude of the white noise. In addition we have 
a set of n coefficients {<f„] describing the power spectrum of the 
GWB, resulting in a total parameter space of {Np{2n + 1) + n) di- 
mensions. 

Those frequencies chosen to be sampled are such that when 
performing the search for the peak in the joint distribution the sam- 
pled set of n coefficients represent the n that contribute the greatest 
total power to the cross correlated spectrum. A more detailed inves- 
tigation into this aspect of the sampling process will be presented 
in an updated version of this paper following the close of the first 
data challenge. 

Table [T] shows a comparison of the values returned for Ai,, a, 
the average amplitude of the white noise cr^^g and the run time re- 
quired for convergence for different values of n, as experienced 
when using a single 16 core Sandy Bridge node on the high perfor- 
mance computer (HPC) 'DARWIN'. We also compare this method 
with that described in vHL2012, running on the same system, 
where we are exploring a 38 dimensional parameter space using the 
MULTINEST algorithm i Feroz, Hobson, & Bridges 2009| |Feroz| 
[& Hobson 2008) to perform parameter estimation (sampling effi- 
ciency set to 0.8 with 2000 live points). 
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Table 1. Comparison of sampling Methods 



Method 


n 


Ah 


a 


sec 


Run Time 


This Paper 


3 


(5.1 ±0.5) X 10-'" 


A.n ± 0.20 


(1.20 ± 0.03) X 10-' 


9 minutes 




4 


(5.2 ± 0.4) X 10-'" 


4.18 ±0.18 


(1.12 ± 0.015) X 10-' 


20 minutes 




5 


(5.09 ± 0.27) X 10-'** 


4.21 ± 0.14 


(1.06 ± 0.03) X 10-' 


35 minutes 




6 


(5.05 ± 0.20) X 10-'" 


4.20 ± 0.12 


(1.040 ± 0.024) X 10-' 


1 hour 


VHL2012 


N/A 


(4.82 ± 0.18) X 10-'" 


4.4 ± 0.08 


(0.97 ± 0.028) X 10-' 


104 hours 



From Table [T] one can see that even with as few as 3 GWB 
power coefficients, the constraints placed on the background are 
only ~ a factor 2 worse than using the method from vHL2012, 
and only 6 coefficients are necessary to match them. However, only 
when using 5+ coefficients does one also constrain the white noise 
to the correct value, with fewer coefficients overestimating its value 
by 10-20%. Run times however are dramatically lower, even in the 
case where 6 coefficients were used, a decrease from ~ 4 days to 
~ 2 hours was observed, representing an improvement of ~ two 
orders of magnitude. 

For the case n = 6 the results of parameterising the white 
noise coeflBcients {Epj can be seen in Fig[T| We find an average 
value for the amplitude of the white noise across all pulsars, iXavg, 
of (1.040 ± 0.024) X 10-', consistent with the data to within 2cr 
confidence levels. 

Fig|2] shows the ID and 2D marginalised posteriors for the 6 
fitted GWB power coefficients {p,j, whilst Fig. [5] shows the power 
spectrum of both the coefficients {a} for each of the 36 pulsars (red 
crosses), and of the GWB power coefficients jpi) (Green crosses) 
with the best fit single power law to the GWB coefficients over- 
laid as a blue dotted line. We find A,, = (5.05 ± 0.2) x 10-'" and 
a = 4.21 ± 0.12, consistent with the injected spectrum at ~ 1 (t 
confidence levels. 

4.2 Open Dataset 2 

The second open data set is designed similarly to the first, with the 
same number of pulsars and observations and the same parameters 
describing the injected GWB signal. The two differences that sep- 
arate the datasets is that the second dataset has uneven sampling in 
the time domain, and the amplitude of the white noise varies per 
pulsar. 

For the case of n = 10 we show the 10 fitted GWB power 
coefficients (green points) with their associated errors, along with 
the best fit power law in Fig|4] We find values for A/, and a of 
(5.4 ± 0.3) X IQ-'" ando' = 4.30 ±0.12, consistent with the injected 
spectrum at ~ 1 cr confidence levels. 



5 CONCLUSIONS 

We have investigated a new method for analysing data from a pulsar 
timing array that results in a speedup of approximately two orders 
of magnitude when compared to methods found in vHL2012 by 
rephrasing the likelihood function to eliminate all matrix-matrix 
multiplications, and costly dense matrix inversions, whilst con- 
straining the parameters of interest within the data to the same level 
so that this new approach scales as 0(log(n)). 

We have shown that in addition to providing a substantial 
speed up, this method offers other advantages; it does not require 



any prior assumptions to be made regarding the shape of the power 
spectrum of the GWB signal within the data. There are also no lim- 
itations placed on the number of parameters to be fit at the point 
of sampling, with reasonable memory usage (~ 1GB) and we have 
shown that it is robust to the precense of irregularly sampled time 
series data, making it suitable for real world applications. 

This is made possible through the use of a 'Guided Hamil- 
tonian Sampler', enabling us to efficiently sample from high- 
dimensional spaces (up to ~ 10'), such that rather than fitting for a 
small number of global parameters that describe a particular model, 
we can simultaneously sample over the set of fourier coefficients 
that describe the spectrum of each pulsar individually, and also the 
coefficients of the GWB itself, along with any additional parame- 
ters of interest, such as the level of the white noise on a per pulsar 
basis. This allows us to make inferences not possible with exist- 
ing methods, and as we will show in a subsequent investigation 
this makes it extremely robust to the presence of bright individual 
sources in the field, or individual pulsars that suffer from extreme 
local effects that might hinder methods of analysis that assume a 
particular model for any cross-correlated signals within the data. 
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Figure 1. Parameterised values for the white noise in each pulsar in open dataset 1 from the IPTA Data Challenge. Each Pulsar has a white noise 
component to their residuals with an amplitude of o",, = lO^^s. Averaging across all pulsars we find an rms value for the white noise of Zavg = -6.983± 
0.010 which is thus consistent with the value in the dataset to within Itr errors. 
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Figure 2. ID and 2D marginalised posteriors for the six GWB Power Spectrum coefficients {p, j. The vertical line in the ID distribution represents the 
power in the injected background at the frequency of that coefficient. Contours in the 2D plots represent 68 and 95 % confidence levels. All coefficients 
are consistent with the injected spectrum at the \ - la confidence level. 
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Figure 3. Log-Log Plot of the parameterised GWB power spectrum in Open Dataset 1 . Tiie red points represent the power at the sampled frequency 
for each of the 36 pulsars, whilst the green points represent the marginalised values of the 6 GWB power coefficients {p, |. The blue dotted line shows 
the best fit power spectrum to the marginalised coefficients, for which A;, = (5.05 ± 0.2) X lO"'** and a = 4.2 ± 0.12, consistent with the injected 
spectrum at ~ 1 cr confidence levels. 
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Figure 4. Log-Log Plot of the parameterised GWB power spectrum in Open Dataset 2. Tiie green points represent the marginalised values of the 10 
GWB power coefficients [pi]. The blue dotted line shows the best fit power spectrum to the marginalised coefficients, for which A;, = (5.4±0.3)x 10"'** 
and a = 4.3 ± 0.12, consistent with the injected spectrum at ~ 1 cr confidence levels. 
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